{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1> <center> Logistic regression </center> </h1>\n",
    "\n",
    "# 1. Model specification\n",
    "\n",
    "Logistic regression corresponding to the following binary classification model\n",
    "$$\n",
    "p(y=1|\\vec{x},\\vec{w})=Ber(y|sigm(\\vec{w}^T\\vec{x})\n",
    "$$\n",
    "where \n",
    "$$\n",
    "sigm(\\eta)=\\frac{1}{1+exp(-\\eta)}\n",
    "$$ \n",
    "is the logistic function.For multi-class logistic regression, the model becomes\n",
    "$$\n",
    "p(y=c|\\vec{x},\\mathbf{W})=\\frac{exp(\\vec{w}_c^T\\vec{x})}{\\sum_{c'=1}^C exp(\\vec{w}_{c'}^T\\vec{x})}\n",
    "$$\n",
    "The negative log-likelihood function of Logistic regression is given as \n",
    "\\begin{align}\n",
    "p(D|\\vec{w}) &=-\\sum_{i=1}^N log \\left[ \\mu_i^{\\mathbb{1}(y_i=1)} \\left(1-\\mu_i)\\right)^{\\mathbb{1}(y_i=0)} \\right] \\\\\n",
    "             &=-\\sum_{i=1}^N \\left[\\mathbb{1}(y_i=1) log \\mu_i+\\mathbb{1}(y_i=0) log(1-\\mu_i) \\right ]\n",
    "\\end{align}\n",
    "where $\\mu_i=sigm(\\vec{w}^T\\vec{x}_i)$. The negative log-likelihood function is also called the **cross-entropy** loss function sometimes.\n",
    "\n",
    "# 2. MLE solution\n",
    "For logistic regression, we can't get the closed form solution anymore. Instead, we need to use an optimization algorithm to compute the MLE solution.\n",
    "The simplest algorithm for unconstrained optimization is **gradient descent**,also known as **steepest descent**. For logistic regression, the negative log-likelihood function is as follows\n",
    "\\begin{align*}\n",
    "NLL(\\vec{w}) &=-\\sum_{i=1}^N \\left[\\mathbb{1}(y_i=1) log \\mu_i+\\mathbb{1}(y_i=0) log(1-\\mu_i) \\right] \\\\\n",
    "\\mu_i &=sigm(\\vec{w}^T\\vec{x_i})\n",
    "\\end{align*}\n",
    "Our task is to minimize the negative log-likelihood function respecting to $\\vec{w}$. We can write the gradient as \n",
    "\\begin{align}\n",
    "\\frac{\\partial NLL(\\vec{w})}{\\partial \\vec{w}} &=-\\sum_{i=1}^N \\lbrace \\mathbb{1}(y_i=1) \\frac{\\partial (log\\mu_i)}{\\partial \\vec{w}} +\\mathbb{1}(y_i=0) \\frac{\\partial (log(1-\\mu_i))}{\\partial \\vec{w}} \\rbrace\\\\\n",
    "                                               &=-\\sum_{i=1}^N \\lbrace \\mathbb{1}(y_i=1) \\frac{1}{\\mu_i}\\frac{\\partial \\mu_i}{\\partial \\vec{w}}+\\mathbb{1}(y_i=0) \\frac{-1}{1-\\mu_i}\\frac{\\partial \\mu_i}{\\partial \\vec{w}}\\rbrace\n",
    "\\end{align}\n",
    "The gradient of the logistic function are given by the following \n",
    "\\begin{align}\n",
    "g &=\\frac{\\partial \\, \\mu}{\\partial \\vec{w}}   \\\\\n",
    "  &=\\frac{\\partial\\, sigm(\\vec{w}^T\\vec{x})}{\\partial \\vec{w}} \\\\\n",
    "  &=\\mu(1-\\mu)\\vec{x}\n",
    "\\end{align}\n",
    "We can get the overall gradient as\n",
    "\\begin{align}\n",
    "\\frac{\\partial NLL(\\vec{w})}{\\partial \\vec{w}} &=-\\sum_{i=1}^N \\lbrace \\mathbb{1}(y_i=1)(1-\\mu_i)\\vec{x}+\\mathbb{1}(y_i=0) (-\\mu_i) \\vec{x} \\rbrace \\\\\n",
    " &=\\sum_{i=1}^N(\\mu_i-y_i)\\vec{x}_i \n",
    "\\end{align}\n",
    "The steepest descent update formular is as follows\n",
    "$$\n",
    "\\vec{w}_{k+1}=\\vec{w}_{k}-\\eta_k \\left( \\frac{\\partial NLL(\\vec{w})}{\\partial \\vec{w}} \\right)_k\n",
    "$$\n",
    "where $\\eta_k$ is the learning rate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The accuracy on validate set of 0 th iteration  is :0.1172 \n",
      "The accuracy on validate set of 1000 th iteration  is :0.8672 \n",
      "The accuracy on validate set of 2000 th iteration  is :0.8594 \n",
      "The accuracy on validate set of 3000 th iteration  is :0.8828 \n",
      "The accuracy on validate set of 4000 th iteration  is :0.8359 \n",
      "The accuracy on validate set of 5000 th iteration  is :0.8438 \n",
      "The accuracy on validate set of 6000 th iteration  is :0.8750 \n",
      "The accuracy on validate set of 7000 th iteration  is :0.8047 \n",
      "The accuracy on validate set of 8000 th iteration  is :0.8125 \n",
      "The accuracy on validate set of 9000 th iteration  is :0.7500 \n"
     ]
    }
   ],
   "source": [
    "import tensorflow as tf\n",
    "import numpy as np\n",
    "from sklearn.model_selection import train_test_split\n",
    "np.random.seed(42)\n",
    "\n",
    "class logistic_regression:\n",
    "    def __init__(self,valid_freq=100):\n",
    "        '''\n",
    "        Initialize the model\n",
    "        \n",
    "        Parameters:\n",
    "            valid_freq:\n",
    "                validate frequence\n",
    "        '''\n",
    "        self.valid_freq=valid_freq\n",
    "    \n",
    "    def fit(self,X_train,y_train,X_valid=None,y_valid=None,max_iters=2000,batch_size=128,learning_rate=0.01):\n",
    "        '''\n",
    "        Fit a logistic model\n",
    "        \n",
    "        Parameters:\n",
    "            X_train:array-like,shape(n_samples,n_features)\n",
    "                train data\n",
    "            y_train:array-like,shape(n_samples,)\n",
    "                train target\n",
    "            X_valid:array-like,shape(n_samples,n_features),optional, default None\n",
    "                validate data\n",
    "            y_valid:array-like,shape(n_samples),optional,default None\n",
    "                validate target\n",
    "            max_iters:int,optional, default 2000\n",
    "                maximum iterations\n",
    "            batch_size:int,optional,default 32\n",
    "                batch size\n",
    "            learning_rate:float,optional,default 0.01\n",
    "                learning rate\n",
    "        '''\n",
    "        with tf.Graph().as_default():\n",
    "            n_features=X_train.shape[1]\n",
    "            n_classes=np.unique(y_train).size\n",
    "            \n",
    "            weight=tf.get_variable(\"weight\",shape=[n_features,n_classes],\n",
    "                                   initializer=tf.random_normal_initializer())\n",
    "            bias=tf.get_variable(\"bias\",shape=[n_classes],\n",
    "                                 initializer=tf.constant_initializer(0.0))\n",
    "            \n",
    "            X_train=X_train.astype(np.float32)\n",
    "            y_train=y_train.astype(np.int64)\n",
    "            \n",
    "            train_dataset=tf.data.Dataset.from_tensor_slices((X_train,y_train))\n",
    "            train_dataset=train_dataset.batch(batch_size).repeat()\n",
    "            \n",
    "            valid_dataset=train_dataset  # valid_dataset must be existed no matter X_valid is None or not\n",
    "            valid_flag= (X_valid is not None and y_valid is not None)\n",
    "            if valid_flag:\n",
    "                X_valid=X_valid.astype(np.float32)\n",
    "                y_valid=y_valid.astype(np.int64)\n",
    "                valid_dataset=tf.data.Dataset.from_tensor_slices((X_valid,y_valid))\n",
    "                valid_dataset=valid_dataset.batch(batch_size).repeat()\n",
    "            \n",
    "            training=tf.placeholder(dtype=tf.bool)\n",
    "            x,y_true=tf.cond(training,lambda:train_dataset.make_one_shot_iterator().get_next(),\n",
    "                            lambda:valid_dataset.make_one_shot_iterator().get_next())\n",
    "            \n",
    "            y_pred=tf.matmul(x,weight)+bias  # y_pred with shape(n_samples,n_classes)\n",
    "            y_most_like=tf.cast(tf.argmax(y_pred,axis=1),tf.int64)\n",
    "            accuracy=tf.reduce_sum(tf.cast(tf.equal(y_true,y_most_like),tf.float32))/ batch_size\n",
    "            loss=tf.nn.sparse_softmax_cross_entropy_with_logits(labels=y_true,logits=y_pred)\n",
    "            train_op=tf.train.GradientDescentOptimizer(learning_rate).minimize(loss)\n",
    "            \n",
    "            init=tf.group(tf.global_variables_initializer(),tf.local_variables_initializer())\n",
    "            with tf.Session() as sess:\n",
    "                sess.run(init)\n",
    "                for i in range(max_iters):\n",
    "                    _=sess.run(train_op,feed_dict={training:True})\n",
    "                    if i%self.valid_freq==0:\n",
    "                        accu=sess.run(accuracy,feed_dict={training:False})\n",
    "                        print(\"The accuracy on validate set of %d th iteration  is :%.4f \" %(i,accu))\n",
    "                        \n",
    "mnist=tf.keras.datasets.mnist\n",
    "(X_train,y_train),(X_valid,y_valid)=mnist.load_data()\n",
    "X_train,X_valid=X_train/255.0,X_valid/55.0\n",
    "X_train=X_train.reshape((X_train.shape[0],-1))\n",
    "X_valid=X_valid.reshape((X_valid.shape[0],-1))\n",
    "\n",
    "model=logistic_regression(valid_freq=1000)\n",
    "model.fit(X_train,y_train,X_valid,y_valid,max_iters=10000,learning_rate=0.01)\n",
    "        "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:tensorflow-cpu-1.11]",
   "language": "python",
   "name": "conda-env-tensorflow-cpu-1.11-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
